1 - Número de SADs preditas não refutadas

Figura 1.1 logito(número de SADs não refutadas) ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.

A figura 1.1 mostra que o padrão geral dos dados não é linear na escala da função de ligação, assim utilizamos GAM. Possibilitamos até 1 intercepto por modelo neutro por sítio de amostragem (MN|SiteCode); para modelos com variáveis de dispersão continua não foi possível um smoother por Sítio de Amostragem e Modelo Neutro

1.1 GLMM binomial

1.1.1 Modelo Cheio e Auditoria

l_md <- vector("list",11)
names(l_md) <- c("k0 1|Site","k0 MN|Site","k0 k0*MN|Site",
                 "kf 1|Site","kf MN|Site",
                 "d 1|Site","d MN|Site","d d*MN|Site",
                 "d/L_plot 1|Site","d/L_plot MN|Site","d/L_plot d/L_plot*MN|Site")
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * k_1.z * MN + I(p.z^2) * k_1.z * MN + (1|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * k_1.z * MN + I(p.z^2) * k_1.z * MN + (MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * k_1.z * MN + I(p.z^2) * k_1.z * MN + (k_1.z*MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * k * MN + I(p.z^2) * k * MN + (1|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[5]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * k * MN + I(p.z^2) * k * MN + (MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[6]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * d.z * MN + I(p.z^2) * d.z * MN + (1|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[7]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * d.z * MN + I(p.z^2)G * d.z * MN + (MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[8]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * d.z * MN + I(p.z^2) * d.z * MN + (d.z*MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[9]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * d_Lplot.z * MN + I(p.z^2) * d_Lplot.z * MN + (1|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[10]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * d_Lplot.z * MN + I(p.z^2) * d_Lplot.z * MN + (MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[11]] <- glmer(cbind(n_nRef,100-n_nRef) ~ p.z * d_Lplot.z * MN + I(p.z^2) * d_Lplot.z * MN + (d_Lplot.z*MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))

Tabela 1.1 AICctab

##                           dAICc   df  weight
## d d*MN|Site                   0.0 22  1     
## d/L_plot d/L_plot*MN|Site   136.6 22  <0.001
## k0 k0*MN|Site              7272.3 22  <0.001
## kf MN|Site                35211.1 123 <0.001
## d MN|Site                 60068.1 15  <0.001
## d/L_plot MN|Site          60179.8 15  <0.001
## k0 MN|Site                61601.1 15  <0.001
## kf 1|Site                 65788.2 121 <0.001
## d 1|Site                  91101.2 13  <0.001
## d/L_plot 1|Site           91273.1 13  <0.001
## k0 1|Site                 93015.1 13  <0.001

O único modelo plausível é aquele com função de ligação logito, a variável d e a estrutura aleatória MN * d |SiteCode.

Figura 1.2 Resíduos Quantílicos do modelo cheio para n_nRef ~ pdMN + p^2dMN + (d*MN|SiteCode)

NOTA: não entendo o motivo de não rodar os resíduos quantílicos quando compilo o Rmarkdown:

mensagem de erro:

Error in if (family$family %in% c(“binomial”,“poisson”,“quasibinomial”, : argumento tem comprimento zero Calls: … withCallingHandlers -> withVisible -> eval -> eval -> simulateResiduals Além disso: warining message: In checkModel(fittedModel) : DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work

1.1.2 Modelo cheio: predições

Figura 1.3 Predito e observado por SiteCode

Observamos que é necessário considerar também um termo quadrático para d. Então sigo com a adição de um termo quadrático para d

Tabela 1.2 R2 condicional e marginal do modelo cheio com termo quadrático para p e d

md_nRef2 <- glmer(cbind(n_nRef,100-n_nRef) ~ 
                     (p.z + I(p.z^2)) * (d.z + I(d.z^2)) * MN + ((d.z + I(d.z^2)) * MN|SiteCode),
                   family = "binomial",data=df_resultados,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
r.squaredGLMM(md_nRef2)

##                   R2m       R2c
## theoretical 0.2153057 0.9548503
## delta       0.2128725 0.9440595

Figura 1.4 Resíduos quantílicos para o modelo cheio com termos quadráticos para p e d. 

Figura 1.5 Predito pelo modelo com termo quadrático para p e d

O modelo cheio com termo quadrático para p e d não foi foi suficiente para descrever o padrão observado nos dados. Hipotetizamos que a razão é a divergência no grau de variação entre os modelos neutros. Assim iremos ajustar modelos para cada conjunto de dados. Para possibilitar a comparação irei manter a estrutura comum das preditoras porem comparando a estrutura aleatória e função de ligação para cada modelo neutro.

1.2 GLMM binomail subset:MN==EE

Modelo cheio

l_md <- vector("list",3)
names(l_md) <- c("1|Site","d|Site","d + d^2|Site")
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (1|SiteCode),
                   family = "binomial",data=filter(df_resultados,MN=="EE"),
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (d.z|SiteCode),
                   family = "binomial",data=filter(df_resultados,MN=="EE"),
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (d.z + I(d.z^2)|SiteCode),
                   family = "binomial",data=filter(df_resultados,MN=="EE"),
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)))

Tabela 1.2.1 AICctab para n_nRef MNEE

##              dAICc   df weight
## d + d^2|Site     0.0 15 1     
## d|Site        2588.4 12 <0.001
## 1|Site       20757.2 10 <0.001

O modelo mais plausível considera um termo linear e quadrático para d na estrutura aletória.

Figura 1.6 Resíduos Quantílicos do glmm cheio para MNEE

Tabela 1.2.2 R2 condicional e marginal modelo cheio nRef

##                     R2m       R2c
## theoretical 0.056151582 0.9244036
## delta       0.006904096 0.1136597

Seleção de Variáveis

# dados
df_md <- filter(df_resultados,MN=="EE")
# modelo global
global_md <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (d.z + I(d.z^2)|SiteCode),
                   family = "binomial",data=df_md,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),na.action = "na.fail")

Tabela 1.2.3 R2 condicional e marginal do modelo global

## Global model call: glmer(formula = cbind(n_nRef, 100 - n_nRef) ~ (p.z + I(p.z^2)) * 
##     (d.z + I(d.z^2)) + (d.z + I(d.z^2) | SiteCode), data = df_md, 
##     family = "binomial", control = glmerControl(optimizer = "bobyqa", 
##         optCtrl = list(maxfun = 1e+05)), na.action = "na.fail")
## ---
## Model selection table 
##     (Int)     d.z   d.z^2    p.z  p.z^2 d.z:p.z d.z:I(p.z^2) I(d.z^2):p.z
## 30  3.856 -0.3820         0.4117 -1.179  0.8548                          
## 32  3.884 -0.3320 0.07890 0.4111 -1.182  0.8581                          
## 62  3.861 -0.4017         0.4087 -1.183  0.8656     0.019820             
## 96  3.885 -0.3384 0.07500 0.3676 -1.184  0.7767                   -0.1387
## 64  3.888 -0.3493 0.07873 0.4084 -1.186  0.8674     0.017350             
## 160 3.889 -0.3312 0.09039 0.4094 -1.187  0.8608                          
## 224 3.955 -0.3349 0.21910 0.3254 -1.252  0.7652                   -0.2256
## 128 3.887 -0.3466 0.07493 0.3664 -1.186  0.7811     0.008296      -0.1385
## 22  2.678 -0.3783         1.1020         0.8500                          
## 192 3.891 -0.3459 0.08725 0.4076 -1.188  0.8671     0.014840             
## 256 3.943 -0.2402 0.24550 0.3323 -1.241  0.7112    -0.094730      -0.2430
## 24  2.706 -0.3298 0.07660 1.1030         0.8531                          
##     d.z^2):I(p.z^2 df    logLik    AICc delta weight
## 30                 11 -6409.908 12841.9  0.00  0.369
## 32                 12 -6409.775 12843.7  1.76  0.153
## 62                 12 -6409.906 12844.0  2.02  0.134
## 96                 13 -6409.357 12844.9  2.95  0.084
## 64                 13 -6409.773 12845.7  3.78  0.056
## 160      -0.011000 13 -6409.774 12845.7  3.78  0.056
## 224      -0.142600 14 -6409.047 12846.3  4.35  0.042
## 128                14 -6409.358 12846.9  4.98  0.031
## 22                 10 -6413.507 12847.1  5.18  0.028
## 192      -0.007982 14 -6409.777 12847.8  5.81  0.020
## 256      -0.168200 15 -6408.992 12848.2  6.27  0.016
## 24                 11 -6413.383 12848.9  6.95  0.011
## Models ranked by AICc(x) 
## Random terms (all models): 
## 'd.z + I(d.z^2) | SiteCode'

Predito para cada Sítio

Para avaliar o ajuste do modelo ao conjunto de dados utilizo funções do pacote MuMIn (REF) que permite calcular a predição do modelo médio para o conjunto original de dados.

Figura 1.7 Predições por sítio de amostragem do modelo médio calculado a partir do modelo global ; pontos: número de SADs neutras refutadas para uma bateria de simulações com um determinada distância média de dispersão; eixo x = distância média de dispersão / largura da área de amostragem; a linha é a probabilidade de não refutar uma SAD neutra segundo a predição média do conjunto de sub-modelos dentro do intervalo de plausibilidade de 7 (Burnham et al 2011)

Predito para novo conjunto de dados

Para avaliar o predito e intervalo de confiança de 95% pelo modelo médio utilizo funções do pacote AICcmodavg (REF).

Figura 1.8 Efeito predito de d e p na Probabilidade de não refutar segundo o modelo médio para MN==EE. No painel superior a predição média, nos paineis inferiores o intervalo de confiança.

1.3 n_nRef subset:MN==EI

Tabela 1.3.1 AICctab n_nRef MNEI

##              dAICc   df weight
## d + d^2|Site     0.0 15 1     
## d|Site        7920.4 12 <0.001
## 1|Site       63552.4 10 <0.001

A função AICcmodavg::modavgPred aceita apenas a função de ligação canonica logito, então não irei comparar outras funções de ligação O modelo mais plausível considera um termo linear e quadrático para d na estrutura aletória. Segue seleção de variáveis.

Tabela 1.3.2 R2 marginal e condicional n_nRef MNEI

##                   R2m       R2c
## theoretical 0.2834906 0.9613372
## delta       0.2808901 0.9525187

Figura 1.9 Resíduos Quantílicos do glmm cheio mais plausível para MNEI

# dados
df_md <- filter(df_resultados,MN=="EI")
# modelo global
global_md <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (d.z + I(d.z^2)|SiteCode),
                   family = "binomial",data=df_md,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),na.action = "na.fail")

Tabela 1.3.3 AICctab n_nRef MNEI delta<7

## Global model call: glmer(formula = cbind(n_nRef, 100 - n_nRef) ~ (p.z + I(p.z^2)) * 
##     (d.z + I(d.z^2)) + (d.z + I(d.z^2) | SiteCode), data = df_md, 
##     family = "binomial", control = glmerControl(optimizer = "bobyqa", 
##         optCtrl = list(maxfun = 1e+05)), na.action = "na.fail")
## ---
## Model selection table 
##     (Int)    d.z  d.z^2       p.z   p.z^2 d.z:p.z d.z:I(p.z^2)
## 79  2.735        -2.147  0.052260 -0.9115                     
## 80  2.935 0.5378 -1.967  0.056100 -0.9083                     
## 208 3.151 0.5308 -2.186 -0.073260 -1.1270                     
## 207 2.952        -2.344 -0.068360 -1.1250                     
## 96  2.923 0.5179 -1.980  0.194700 -0.9040  0.3705             
## 112 3.154 0.8376 -1.981 -0.009226 -1.1350              -0.3201
## 224 3.151 0.5320 -2.184  0.063220 -1.1270  0.3650             
## 240 3.214 0.6970 -2.134 -0.074240 -1.1910              -0.1693
## 128 3.070 0.7209 -1.980  0.109200 -1.0510  0.2517      -0.2028
## 256 3.109 0.4198 -2.219  0.088050 -1.0840  0.4312       0.1148
##     I(d.z^2):p.z d.z^2):I(p.z^2 df    logLik    AICc delta weight
## 79       -1.1980                11 -6288.847 12599.8  0.00  0.180
## 80       -1.1980                12 -6287.921 12600.0  0.17  0.165
## 208      -1.0710         0.2145 13 -6287.193 12600.6  0.74  0.124
## 207      -1.0830         0.2025 12 -6288.288 12600.7  0.91  0.114
## 96       -1.0740                13 -6287.311 12600.8  0.98  0.110
## 112      -1.1340                13 -6287.345 12600.9  1.04  0.107
## 224      -0.9524         0.2140 14 -6286.677 12601.6  1.74  0.075
## 240      -1.0710         0.1598 14 -6287.097 12602.4  2.58  0.050
## 128      -1.0740                14 -6287.169 12602.5  2.72  0.046
## 256      -0.9315         0.2507 15 -6286.649 12603.5  3.71  0.028
## Models ranked by AICc(x) 
## Random terms (all models): 
## 'd.z + I(d.z^2) | SiteCode'

Predito para cada Sítio

Para avaliar o ajuste do modelo ao conjunto de dados utilizo funções do pacote MuMIn (REF) que permite calcular a predição do modelo médio para o conjunto original de dados.

Figura 1.10 Predito por SiteCode a partir do modelo médio para MNEI.

Figura 1.11 Efeito predito de d e p na Probabilidade de não refutar segundo o modelo médio para MN==EE. No painel superior a predição média, nos paineis inferiores o intervalo de confiança.

1.4 Figura Final

Figura 1.12 Probabilidade de não refutar uma SAD neutra em função do modelo neutro (EE e EI), proporção de cobertura vegetal (p) e distância média de dispersão (d)

2 diff_S = (S_obs - S_MN)/S_obs

Padrão Geral

Figura 2.1 diff_S (y-axis), p (x-axis).

2.1 MNEE

l_md <- vector("list",9)
names(l_md) <- c("1-k 1|Site","1-k I(1-k)|Site",
                 "k 1|Site",
                 "d 1|Site","d d|Site",
                 "d/L 1|Site","d/L I(d/L)|Site",
                 "1 1|Site",
                 "1")
df_md <- df_resultados %>% filter(MN=="EE")
l_md[[1]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * k_1.z + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[2]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * k_1.z + (k_1.z|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[3]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * k + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[4]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * d.z + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[5]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * d.z + (d.z|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[6]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * d_Lplot.z + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[7]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * d_Lplot.z + (d_Lplot.z|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[8]] <- lmer(diff_S0 ~ 1 + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[9]] <- lm(diff_S0 ~ 1,
                  data=df_md,na.action = "na.fail")
AICctab(l_md,weights=T)

##                 dAICc df weight
## 1                 0.0 2  0.928 
## 1 1|Site          5.1 3  0.072 
## 1-k 1|Site       73.3 8  <0.001
## 1-k I(1-k)|Site  73.6 10 <0.001
## d/L 1|Site       73.9 8  <0.001
## d 1|Site         74.0 8  <0.001
## d/L I(d/L)|Site  75.1 10 <0.001
## d d|Site         75.6 10 <0.001
## k 1|Site        711.9 62 <0.001

Alguns modelos utilizados para estimar diff_S apresentaram “singularidade” (estimativas não diferem de zero), dentre eles o único modelo plausível. Isso indica que modelo com a estrutura aleatória mais simples devem ser ajustados aos dados, assim vou descartar aqueles com inclinação para a variável de dispersão por Sítio.

l_md <- vector("list",6)
names(l_md) <- c("1-k 1|Site",
                 "k 1|Site",
                 "d 1|Site",
                 "d/L 1|Site",
                 "1 1|Site",
                 "1")
df_md <- df_resultados %>% filter(MN=="EE")
l_md[[1]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * k_1.z + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[2]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * k + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[3]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * d.z + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[4]] <- lmer(diff_S0 ~ (p.z + I(p.z^2)) * d_Lplot.z + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[5]] <- lmer(diff_S0 ~ 1 + (1|SiteCode),
                   data=df_md,na.action = "na.fail")
l_md[[6]] <- lm(diff_S0 ~ 1,
                   data=df_md,na.action = "na.fail")
AICctab(l_md,weights=T)

##            dAICc df weight
## 1            0.0 2  0.928 
## 1 1|Site     5.1 3  0.072 
## 1-k 1|Site  73.3 8  <0.001
## d/L 1|Site  73.9 8  <0.001
## d 1|Site    74.0 8  <0.001
## k 1|Site   711.9 62 <0.001

Nenhum destes modelos apresentou singularidade, o modelo nulo sem estrutura aleatória foi o único plausível. Segue gráfico diagnóstico:

Figura 2.2 Resíduos Quantílicos do modelo mais plausível para diff_S subset=MNEE (diff_S ~ 1)

Tabela 2.3 Sumário do modelo mais plausível para diff_S MNEE

## 
## Call:
## lm(formula = diff_S0 ~ 1, data = df_md, na.action = "na.fail")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.080599 -0.009762  0.000602  0.010029  0.054253 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.003539   0.000342  -10.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01552 on 2059 degrees of freedom

MNEI

Figura 2.3 diffS_EI ~p*k

Vou iniciar apenas com as preditoras de interesse (p e variáveis de dispersão). Devido a clara ausência de simetria de diff_S desloquei a distribuição apenas para valores positivos (figura 2.3 primeiro quadro) e utilizarei distribuição Gamma para ajustar os modelos com função de ligação ‘log’.

Figura 2.4 diff_S por sítio de amostragem e d

Seleção do Modelo cheio: GLMM

df_md <- df_resultados %>% filter(MN=="EI")
# df_md$diff_S %>% summary
l_md <- vector("list",6)
names(l_md) <- c("1|Site","d|Site",
                 "p2*k log",
                 "p2*d2 id 1|Site",
                 "p2*d2 id d|Site",
                 "p2*d2 id d2|Site")
# modelos mínimos
l_md[[1]] <- glmer(diff_S ~ p.z*d.z + (1|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
l_md[[2]] <- glmer(diff_S ~ p.z*d.z  + (d.z|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
# com p2
## erro
l_md[[3]] <- glmer(diff_S ~ (p+I(p.z^2))*k + (1|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
v_vars <- getME(l_md[[3]],c("theta","fixef"))
l_md[[3]] <- update(l_md[[3]],start=v_vars,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
l_md[[3]] <- glmer(diff_S ~ I(p.z^2)*k + (1|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
# com p2 e d2
#######################
# Model failed to converge with max|grad| = 0.00556207 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
#  - Rescale variables?
l_md[[4]] <- glmer(diff_S ~ (p.z + I(p.z^2))*(d.z + I(d.z^2)) + (1|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
v_vars <- getME(l_md[[4]],c("theta","fixef"))
l_md[[4]] <- update(l_md[[4]],start=v_vars,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#######################
#######################
# Model failed to converge with max|grad| = 0.0608783 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
#  - Rescale variables?
l_md[[5]] <- glmer(diff_S ~ (p.z + I(p.z^2))*(d.z + I(d.z^2)) + (d.z|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
v_vars <- getME(l_md[[5]],c("theta","fixef"))
l_md[[5]] <- update(l_md[[5]],start=v_vars,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#######################
# Model failed to converge with max|grad| = 0.336616 (tol = 0.002, component 1)Model is nearly unidentifiable: very large eigenvalue
#  - Rescale variables?
l_md[[6]] <- glmer(diff_S ~ (p.z + I(p.z^2))*(d.z + I(d.z^2)) + (d.z + I(d.z^2)|SiteCode),
                   data=df_md,family="Gamma"(link="log"))
v_vars <- getME(l_md[[6]],c("theta","fixef"))
l_md[[6]] <- update(l_md[[6]],start=v_vars,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5)))
#######################

##                  dAICc  df weight
## p2*k log            0.0 42 1     
## p2*d2 id d2|Site 5496.0 16 <0.001
## p2*d2 id d|Site  5810.5 13 <0.001
## p2*d2 id 1|Site  5842.6 11 <0.001
## d|Site           8272.9 8  <0.001
## 1|Site           8300.0 6  <0.001

Não foi possível ajustar os modelos com estrutura fixa mais complexa. Vou avaliar se GAMM oferecem melhor ajuste aos dados

GAMM cheio

l_md <- vector("list",4)
names(l_md) <- c("tp by=k","cr by=k",
                 "tp ti","cr ti")
l_md[[1]] <- gam(diff_S ~ s(p.z,by=k,bs="tp") + 
                   s(SiteCode,bs="re"),
                 data=df_md,family = "Gamma"(link="log"))
l_md[[2]] <- gam(diff_S ~ s(p.z,by=k,bs="cr") + 
                   s(SiteCode,bs="re"),
                 data=df_md,family = "Gamma"(link="log"))
l_md[[3]] <- gam(diff_S ~ s(p.z,bs="tp") + s(d.z,bs="tp") + ti(p.z,d.z,bs=c("tp","tp")) +
                   s(SiteCode,bs="re"),
                 data=df_md,family = "Gamma"(link="log"))
l_md[[4]] <- gam(diff_S ~ s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                   s(SiteCode,bs="re"),
                 data=df_md,family = "Gamma"(link="log"))

##         dAICc   df               weight
## cr ti       0.0 127.257591150728 1     
## tp ti    1081.0 126.582457404983 <0.001
## cr by=k 13007.1 115.96208485403  <0.001
## tp by=k 13007.1 115.962163752688 <0.001

O único modelo cheio plausível considera as variáveis p e d (um smoother para cada mais um tensor para ambas) e 1 intercepto por Sítio de Amostragem.

Figura 2.5 Gráficos Diagnostico do gamm cheio mais plausível para diff_S MN==EI

De maneira geral o modelo esta fazendo um bom ajuste, contudo alguns sítios distoam do padrão geral dos dados. Há duas alternativas: i) possibilitar um smoother por sítio de amostragem com um parâmetro de penalização comum; ou ii) remover sítios outliers.

l_md <- vector("list",2)
names(l_md) <- c("1|Site","d|Site")
l_md[[1]] <- gam(diff_S ~ s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                   s(SiteCode,bs="re"),
                 data=df_md,family = "Gamma"(link="log"))
l_md[[2]] <- gam(diff_S ~ s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                   s(d.z,SiteCode,bs="fs",xt="cr"),
                 data=df_md,family = "Gamma"(link="log"))

##        dAICc df              
## d|Site   0.0 472.606086811532
## 1|Site 858.6 127.257591150728

Figura 2.6 Gráfico Diagnóstico gamm cheio com segunda estrutra aleatória diff_S MM==EI

Figura 2.7 gratia::draw() gamm mais plausível com estrutura similar a (d|Site)

Este modelo apresenta problema de concurvidade. Assim vou refazer a comparação com configuração adequada:

df_md <- filter(df_resultados,MN=="EI")
md_gam.diff_S.EI__update <- gam(diff_S ~ s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                   s(d.z,SiteCode,bs="fs",xt=list(bs="cp"), m=1),
                 data=df_md,family = "Gamma"(link="log"))

Figura 2.8 Update do modelo da figura 2.7 com o smoother da estrutura aleatória penalizada pela 1a derivada para diminuir a concurvidade (Perdersen et al. 2019)

O modelo diminuiu a concurvidade da estimativa.

Figura 2.9 Appraise do gamm cheio atualizado (smoother de estrutura aleatória: m=1).

Remoção de Outliers

Remoção sítios com observações fora do intervalo [-1.5;1.5] dos resíduos deviance.

# dados completos
df_md <- df_resultados %>% filter(MN=="EI")
df_md$resDev <- residuals.gam(md_gam.diff_S.EI__update)
# sítios considerados outliers
v_SiteOut <- df_md %>% filter(resDev < -0.15 | resDev > 0.15) %>% .$SiteCode %>% unique
# update da comparação
l_md <- vector("list",2)
names(l_md) <- c("1|Site","d|Site")
l_md[[1]] <- gam(diff_S ~ s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                   s(SiteCode,bs="re"),
                 data=filter(df_md,!(SiteCode %in% v_SiteOut)),
                 family = "Gamma"(link="log"))
l_md[[2]] <- gam(diff_S ~ s(p.z,bs="cr") + s(d.z,bs="cr") + ti(p.z,d.z,bs=c("cr","cr")) +
                   s(d.z,SiteCode,bs="fs",xt=list(bs="cp"),m=1),
                 data=filter(df_md,!(SiteCode %in% v_SiteOut)),
                 family = "Gamma"(link="log"))

##        dAICc df               weight
## d|Site   0.0 324.888880261917 1     
## 1|Site 472.9 123.947848190999 <0.001

Figura 2.10 Gráficos diagnóstico do gamm cheio mais plausível sem outliers

A remoção de outliers não melhorou o ajuste do modelo.

Figura 2.11 Efeitos estimados para o gamm cheio diffS_EI sOut

Seleção de Variáveis

Figura Final

3 U

Padrões Gerais

figura 4.1 U ~ padrões gerais

Figura 4.2 U ~ d (~SiteCode)

GAMM

Modelo Cheio: comparação estrutura aleatória

Variável de dispersão: d (como para n_nRef)

1. Na escala log a variável parece ser simetrica, então vou comparar a distribuição normal e gamma;
   
2. É possível apresentar 3 estruturas: a) 1|Site; b) d|Site, penalizacao comum; c) d|Site, penalizacao por Site [que não vou incluir a priori pelo custo computacional]
   
3. pelo menos 2 smoother type (tp,cr)

Figura 4.3 Graf Diag Gamm cheio mais plausível

4 Comparação n_nRef e diff_S

Figura 3.1 Observado (pontos) e predito (linhas) para diffS e n_nRef (segunda linha no facet) em função de d*MN, por sítio de amostragem (primeira linha no facet).

Figura 4.2 |diffS| ~ n_nRef * k, subset=MNEI

Figura 3.3 n_nRef ~ |diffS| (~SiteCode), subset=MNEI

Comparação funções para o teste KS

df_testeKS <- df_auditoria %>% select(SAD_obs.name,SAD_MN.name,MN,k,rep,ordem,refID,SiteCode,txt.name,S_obs,p,d) %>% 
  group_by(SAD_obs.name) %>% nest
##### Resultados por Réplica #####
df_testeKS$resultados <- vector("list",length = nrow(df_testeKS))
for(i in 1:nrow(df_testeKS)){
  # i <- 1
  df_ <- df_testeKS[i,]
  v_SAD.obs <- read.csv(df_$SAD_obs.name,header = TRUE,as.is = TRUE) %>% filter(species.correct != "Mortas") %>% 
    .$N %>% sort()
  df_predicao <- as.data.frame(df_$data[[1]])
  f_KSeS <- function(v_OBS = v_SAD.obs, path_SAD.MN){
    # path_SAD.MN <- df_predicao[1,"SAD_MN.name"]
    v_SAD.predita <- as.integer(read.csv(path_SAD.MN,header = TRUE,as.is = TRUE)$SAD_predita)
    a <- disc_ks_test(x=v_SAD.predita,
                      y=ecdf(v_SAD.obs),exact = TRUE)
    a <- data.frame(KS.D = a$statistic, KS.p = a$p.value)
    a$S_SAD.predita <- length(v_SAD.predita)
    a$S_SAD.obs <- length(v_SAD.obs) 
    return(a)
  }
  # f_KSeS(path_SAD.MN = df_predicao$SAD_MN.name[3000])
    l_ <- vector("list",length = nrow(df_predicao))
  for(j in 1:nrow(df_predicao)){
    l_[[j]] <- f_KSeS(path_SAD.MN = df_predicao$SAD_MN.name[j])
  }
  df_testeKS$resultados[[i]] <- rbind.fill(l_)
}
df_replicas <- df_testeKS %>% unnest(cols = c(data,resultados)) %>% as.data.frame() %>% 
  rename(D_KSgeneral = KS.D, p.valor_KSgeneral=KS.p)
df_replicas <- inner_join(x=df_replicas,
                          y=select(df_auditoria,
                                   SAD_obs.name,k,rep,SiteCode,KS.D,KS.p),
                          by=c("SAD_obs.name","k","rep","SiteCode"))
write.csv(df_replicas,file = "~/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/simulacao/resultados/df_replicas__KSgeneral.csv",row.names = FALSE)

Visualização das SADs preditas e observadas

Inspeção individual das SADs

Estratégia

1. estrutura do tipo:

cbind(n_nRef,100-n_nRef) ~ log(modulo_diffS) * d + ( log(modulo_diffS) * d | SiteCode)

2. expectativa:

1. esperavamos que se o teste é adequado então deve existir um erro de estimativa da riqueza máximo no qual a probabilidade de não se refutar uma SAD se torna muito diminuta.
   
2. e uma vez que não é suficiente apresentar boa estimativa da riqueza para apresentar boa congruência, então esperava que a nuvem de pontos pudesse ser aproximada por um triangulo retangulo com angulo reto na origem.

##                                    dAICc   df weight
## log(modulo_diffS) * d.z | SiteCode     0.0 27 1     
## log(modulo_diffS) + d.z | SiteCode   794.7 16 <0.001
## d.z | SiteCode                      1727.8 12 <0.001
## log(modulo_diffS) | SiteCode       16432.1 9  <0.001
## 1 | SiteCode                       40804.1 7  <0.001